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FOREWORD 

This  report  describes  the  implementation  of  the  Mie-Gmeneisen  and  P-a  equations  of 
state  in  the  DYSMAS  code.  The  motivation  for  these  enhancements  is  to  provide  an  equation 
of  state  that  treats  partially  saturated  sand.  The  pores  in  the  P-a  model  are  used  to  simulate 
the  sand  pore  volume  free  from  water,  while  the  Mie-Grueneisen  equation  of  state  accounts 
for  the  water-grain  mixture.  Additionally,  the  Mie-Grueneisen  equations  of  state  provides  a 
means  of  modeling  the  equations  of  state  of  metals. 
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programs. 


Approved  by: 


MARY  E  LACEY,  Acting  Head 
Systems  Research  &  Technology  Department 


NSWCDD/rR-95/107 


CONTENTS 

Chapter  Page 

1  INTRODUCTION .  1 

2  MIE-GRUENEISEN  EQUATION  OF  STATE .  2 

EQUATION .  2 

SOUND  SPEED .  3 

CHANGE  OF  STATE .  4 

DELTA  FORMULATION .  5 

DYSMAS  IMPLEMENTATION .  6 

3  P-a  EQUATION  OF  STATE .  7 

DESCRIPTION .  7 

ELASTIC  AND  PLASTIC  a  RELATIONS .  7 

COMPUTING  p  AND  p .  9 

CHANGES  OF  STATE .  9 

DYSMAS  IMPLEMENTATION .  10 

4  NUMERICAL  RESULTS .  12 

VALIDATION  TESTS .  12 

APPLICATION  OF  THE  P-a  TO  SAND  TESTS .  14 

5  SUMMARY  AND  CONCLUSIONS .  16 

REFERENCES .  17 

GLOSSARY .  18 

DISTRIBUTION .  (1) 

Appendix 

A  MIE-GRUNEISEN  EQUATION  OF  STATE  ROUTINES .  A-1 

B  P-a  EQUATION  OF  STATE  ROUTINES .  B-1 

C  USER  MANUAL  PAGES .  C-1 

D  STORAGE  ARRAYS  FOR  EQUATION  OF  STATE  PARAMETERS  ...  D-1 

E  COMPUTING  THE  P-a  PARAMETERS .  E-1 


V 


NSWCDD/TR-95/107 


ILLUSTRATIONS 

Figure  Page 

1  P-a  ELASTIC  PLASTIC  CURVES .  20 

2  POSSIBLE  P-a  SOLUTIONS  of  p(a,  p,  e)  =  P{a^in,p) .  20 

3  NON-LINEAR  P-a  ELASTIC-PLASTIC  CURVES  FOR  A  SAND  LIKE 

MATERIAL .  21 

4  COMPUTED  AND  EXACT  DENSITY  FOR  THE  MIE-GRUNEISEN 

RIEMANN  PROBLEM .  22 

5  COMPUTED  AND  EXACT  DENSITY  FOR  THE  P-a  RIEMANN 

PROBLEM .  22 

6  PRESSURE  RESPONSE  OF  FULLY  SATURATED  SAND  TO  A 

GAMMA  LAW  EXPLOSION  BUBBLE .  23 

7  VELOCITY  RESPONSE  OF  FULLY  SATURATED  SAND  TO  A 

GAMMA  LAW  EXPLOSION  BUBBLE .  24 

8  PRESSURE  RESPONSE  OF  PARTIALLY  SATURATED  SAND 

(ao  =1.052808)  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE .  25 

9  VELOCITY  RESPONSE  OF  PARTIALLY  SATURATED  SAND 

(ao  =1.052808)  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE .  26 

10  COMPUTE  a  FOR  PARTIALLY  SATURATED  SAND  (ao  =1.052808) 

TO  A  GAMMA  LAW  EXPLOSION  BUBBLE .  27 

11  PRESSURE  HISTORIES  AT  THREE  POINTS  FOR  PARTIALLY 
SATURATED  SAND  (ao  =1.052808)  IN  RESPONSE  TO  A  GAMMA 

LAW  EXPLOSION  BUBBLE .  28 

12  VELOCITY  HISTORIES  AT  THREE  POINTS  FOR  PARTIALLY 
SATURATED  SAND  (ao  =1.052808)  IN  RESPONSE  TO  A  GAMMA 

LAW  EXPLOSION  BUBBLE .  29 

13  INFLUENCE  OF  FCT  PARAMETERS  ON  THE  COMPUTED 

PRESSURE  FOR  A  GAMMA  LAW  EXPLOSION  BUBBLE .  30 

14  INFLUENCE  OF  FCT  PARAMETERS  ON  THE  COMPUTED 

VELOCITIES  FOR  A  GAMMA  LAW  EXPLOSION  BUBBLE .  31 

15  GRAY  SCALE  PRESSURE  PLOT  FOR  A  TNT  EXPLOSION  AT  THE 

WATER-SAND  INTERFACE .  32 

16  GRAY  SCALE  a  PLOT  FOR  AN  UNDERWATER  TNT  EXPLOSION 

AT  THE  WATER-SAND  INTERFACE .  33 


VI 


NSWCDDyTR-95/107 


ILLUSTRATIONS  (CONTINUED) 

Figure  Page 

17  PRESSURE  ALONG  THE  CENTERLINE  OF  THE  EXPLOSION  OF 

FIGURES  15  AND  16 .  34 

18  PRESSURE  ALONG  THE  z  =  -550  LINE  OF  THE  EXPLOSION  OF 

FIGURES  15  AND  16 .  35 

19  PRESSURE  ALONG  THE  z  =  -450  LINE  OF  THE  EXPLOSION  OF 

FIGURES  15  AND  16 .  36 

20  COMPUTED  AND  MEASURED  VELOCITIES  FOR  THE  100% 

SATURATED  SRI  TEST .  37 

21  COMPUTED  AND  MEASURED  VELOCITIES  FOR  THE  95% 

SATURATED  SRI  TEST .  38 

22  COMPUTED  AND  MEASURED  VELOCITIES  FOR  THE  78% 

SATURATED  SRI  TEST .  39 

23  SAMSI  TEST  SETUP .  40 

24  COMPUTED  PRESSURES  FOR  THE  SAMSI  TEST  AT  734  fiSECs  ...  41 

25  COMPUTED  PRESSURES  FOR  THE  SAMSI  TEST  AT  789  ^SECs  ...  41 

26  COMPUTED  PRESSURES  FOR  THE  SAMSI  TEST  AT  844  //SECs  ...  42 

27  COMPUTED  PRESSURES  FOR  THE  SAMSI  TEST  AT  953  fxSECs  ...  42 

28  COMPUTED  AND  MEASURED  SAND  PRESSURES  FOR  THE 

SAMSI  TEST .  43 

29  COMPUTED  CENTERLINE  DISPLACEMENT  HISTORY  FOR  THE 

SAMSI  TEST .  44 

30  COMPURED  AND  MEASURED  TARGET  DEFORMATION  FOR  THE 

SAMSI  TEST .  45 

TABLES 

Table  Page 

1  MIE-GRUENEISEN  AND  P-a  PARAMETERS  FOR  THE 

TEST  CASES .  12 


Vll 


NSWCDD/rR-95/107 


CHAPTER  1 
INTRODUCTION 

i  The  P-a  equation  of  state  was  original  proposed  by  Herrman^  to  model  porous  material. 

This  equation  of  state  assumes  that  a  porous  material  wUl  initially  behave  elastically  when 
I  a  stress  is  applied.  As  the  stress  is  increased,  the  pores  in  the  material  start  to  crushed, 

which  is  an  irreversible  phenomenon  and  leads  to  plastic  compression.  Unloading  of  the 
partially  crushed  porous  material  follows  a  new,  elastic  curve  that  depends  on  the  maximum 
stress  achieved  by  the  material.  At  sufficiently  high  stress,  all  pores  are  eliminated,  and  the 
equation  of  state  behavior  is  that  of  the  original  solid  material.  This  P-a  model  has  been  used 
in  its  original  form  in  the  WONDY  code^  and  a  modified  version  of  it  has  been  included 
'  in  the  CTH  code^.  This  modification  was  proposed  by  Carroll  and  Holt"^  and  is  designed  to 

improve  the  consistency  of  the  P-a  model. 

Although  the  purpose  of  the  P-a  equation  of  state  is  to  treat  porous  material,  the  focus  in 
this  report  is  sand.  Sand  consists  of  sand  grains,  water  and  air  or  void.  This  last  component 
will  be  absent  if  the  sand  is  fully  saturated.  However,  in  most  natural  situations  this  will  not 
be  the  case.  The  solid  part  of  the  P-a  model  represents  the  sand  grain-water  mixture,  while 
the  porous  part  accounts  for  the  air-void  content. 

It  is  well  known  that  sand  supports  a  shear  stress  and  several  Effective  Stress  Models 
have  been  formulated  to  represent  this  characteristic.^  Although  the  P-a  equation  of  state 
cannot  directly  treat  this  phenomenon,  the  shear  strength  of  sand  can  be  addressed  by  the 
application  of  a  deviatoric  stress  model  in  conjunction  with  the  P-a  equation  of  state. 

The  Mie-Grueneisen  equation  of  state  is  used  to  model  the  solid  state  associated  with 
the  P-a  model.  Chapter  2  outlines  the  Mie-Grueneisen  model  and  its  incorporation  into 
DYSMAS.  Chapter  3  discusses  the  P-a  equation  of  state  and  its  implementation  while  Chapter 
4  presents  numerical  examples.  A  summary  of  the  report  is  provided  in  Chapter  5. 
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CHAPTER  2 

MIE-GRUENEISEN  EQUATION  OF  STATE 


The  form  of  the  Mie-Grueneisen  equation  of  state  described  in  Reference  3  has  been 
employed. 


EQUATION 

The  Mie-Grueneisen  equation  of  state  provides  p  as  a  function  of  p  and  e.  The  basic 
form  is; 


p  =  Pr(/>)  +  ro/5o(e-er(/9))  (1) 

where  F  =  and  Tp  =  Fop©.  The  subscript  r  refers  to  the  reference  state  which  is  an 
isentrope  passing  through  the  reference  point  (po,  po,  Co)-  The  reference  state  is  defined  by: 


,  Po  ,  cIp^Y 

Cr  =  eo  +  P  +0/^1  c  ^ 
po  2(1  -  Sp) 

(2) 

Poclp  (,,  dY  Y  \ 

2(1  -  S^)  +  (1  -  5//)  j  ■ 

(3) 

Here: 

(4) 

P 

and 

5 

(5) 

Jfc=0 


The  ak  coefficients  are  defined  by: 

S  1 

ao  =  1;  cn  =  ajfc  =  7— +  kS)a],-i  -  Fo5aA:_2]  (6) 

O  K  Z 

where  S  is  the  coefficient  relating  shock  velocity  and  particle  velocity  and  Fo  is  the  Grueneisen 
parameter  at  the  initial  density  po-  Note  the  definition  of  pr  follows  from  er  via: 

Pt  _  //y. 


since  the  reference  state  is  isentropic. 
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SOUND  SPEED 


An  additional  variable  required  by  DYSMAS  is  sound  speed.  This  is  derived  as  follows: 


2  _  f  MPi  e) 
V  dp 


c  = 


2  _  ( dp\  ( dp\  f  dp\  / de' 

~  \dp)Xdp) \de) 

^  =  f  L. 

\dp)e  \d^) pP^  ' 


(8) 


Using  Equations  (1)  and  (7)  : 


dp 

de 


—  ^OpO 


( dp\  _ 

\Pp)~ 


dpr  ^oPo 


e  dp  p- 

Substituting  these  definitions  into  Equation  (8)  yields  the  expression  for  sound  speed: 

2  dpr  TopOf  \ 

c  =  ^  +  —rip  -  Pr)  ■ 
dp  P^ 


(9) 

(10) 


(11) 


To  evaluate  it  is  necessary  to  compute  dprjdp.  Differentiating  Equation  (3)  produces: 


Plcj 


dpr  Po  dpr 

dp  p  dp  2(1  —  Sp)p^ 


2Y 


2p{2  -  Sp)dY  2 

I  /t  /nr  \  7  I 


(1  -  SpY  (1  -  Sp)  dp 


<fiY 

dp^ 


(12) 


Rearranging  Equation  (11)  and  substituting  Equation  (10)  yields  an  additional  relation 
which  is  applied  in  the  DYSMAS  equation  of  state  routines  is: 


/dp\ 

2  ^oPoP 

II 

^  9 

P  J 

(13) 
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CHANGES  OF  STATE 

The  DYSMAS  code  requires  isentropic  and  Hugoniot  changes  of  state  calculations:  given 
an  initial  state  ipi,  ei)  and  a  second  density,  p2,  compute  the  corresponding  energy  e2  which 
is  on  the  original  isentrope  or  Hugoniot. 

For  an  isentropic  process 


de  = 


pdp 


^  -  e,) 


dp 


Since  the  reference  state  is  isentropic, 


Equation  (14)  reduces  to: 


der  =  ^dp 
P 


d{e  -  er)  _  -p  dp 
(e-er) 


Integrating  this  equation  produces  the  isentropic  change  of  state  formula. 

'^oPo{p2  —  Pi) 


e2  =  er{p2)  +  [ei  -  er{pi)]exp(^ 


P2pl 


(14) 


(15) 


(16) 


(17) 


The  Hugoniot  change  of  state  relation  follows  directly  from  the  Hugoniot  condition 

2/3iP2(e2  -  ei) 


{P2-P1)  • 

Eliminating  e  from  this  equation  via  Equation  (1)  and  solving  for  p2  produces 

,  P1P2C  -  Pl{p2  -  Pi) 

P2  =  W  +  — - . 

2  Topo 


(18) 


where: 


(19) 


C  =  <  er(/?2)  “  ^r{p\)  + 


[Pt{P2)  -  Pt{Pi)] 

^OpO 
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DELTA  FORMULATION 


For  compatibility  with  the  DYSMAS  code,  Equation  (1)  is  written  in  delta  form  which 
references  pressure,  density,  and  energy  to  an  ambient  state  {poo,Poo,^oo)-  Each  material 
has  its  own  density  and  energy  ambient  state,  while  a  common  ambient  pressure  is  used  for 
all  materials.  This  leads  to  the  following  definitions  for  the  delta  variables: 

Ap  =  p- Poo 

Apr  =  Pr  -  Poo 

Ae  =  e  -  too  (21) 

/^tf  tj'  ““  too 

Ap  =  p-  poo 


Recasting  fx  in  delta  form  produces 


^P  ^  {Poo  -  Po) 
P  P 


(22) 


The  delta  form  diminishes  the  truncation  error  in  expressions  which  contain  the  difference 
between  to  state  variables.  This  issue  is  of  concern  since  DYSMAS  is  a  single  precision  code. 
The  equations  necessary  to  implement  the  Mie-Grueneisen  equation  of  state  in  DYSMAS 
are  re-cast  as  follows: 


Eq.  (1)  :  Ap  =  Apr  -h  FoPoIAe  -  Ae^] , 


Eq.{ll):  =  ^+^^{Ap-Apr)  , 


dp 


(i?), 


2  Vopo{,Apr  -f  Poo) 


Eq.  (17)  :  Ae,  =  Ae,(p2)  +  [Ae,  -  Ae.(/,i)]expf  _ 

V  PlP2  J 

Eq.  (19)  :  8p=  -  Api)^  C  =  8er-\-  — ^  ^ 


Sp  _  p\p2 

2  tX 


TqPo 


(23) 

(24) 

(25) 

(26) 

(27) 


5 


NSWCDD/rR-95/107 


DYSMAS  IMPLEMENTATION 

The  DYSMAS  code  requires  equation  of  state  routines  which  perform  the  6  different 
functions  listed  below: 

•  Function  fi:  Compute  Apijcf  given  A/)i,Aei, 

•  Function  £2:  Compute  Ap2,Ae2,c|  given  Ajoi,Aei,A^2  and  assuming  an  isentropic 
change  of  state  between  Api  and  Ap2, 

•  Function  £3:  Compute  Ap2,Ae2,  cj  given  Api,  Aei,  Ap2  and  assuming  a  change  of  state 
along  the  Hugoniot  between  Api  and  Ap2, 

•  Function  £4:  Compute  Api,cl  given  Api,  Aei, 

•  Function  fs:  Compute  Ap2,Ae2,c^  given  Api,Aei,Ap2  and  assuming  an  isentropic 
change  of  state  between  Api  and  Ap2, 

•  Function  fe:  Compute  Ap2,  Ae2,  c|  given  Api,  Aei,  Ap2  and  assuming  a  change  of  state 
along  the  Hugoniot  state  between  Api  and  Ap2. 

A  description  of  these  subroutines  is  provided  in  Appendix  A. 
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CHAPTER  3 

P-a  EQUATION  OF  STATE 


DESCRIPTION 


The  basic  form  of  the  P-a  equation  is: 


V  =  p{ot,p,e) 


a 


(28) 


where 

a  =  ^  .  (29) 

P 

Here  the  subscript  s  denotes  properties  of  the  underlying  solid  material  and  the  variables 
that  are  not  subscripted  pertain  to  the  porous  media.  The  condition  of  the  porous  media  is 
dictated  by  the  function 


Ol-A{p,pt,Pmax)  (30) 

which  is  sketched  in  Figure  1  and  has  multiple  elastic  branches.  Here  pt  is  the  material 
derivative  of  p;  positive  values  indicate  loading,  while  negative  ones  indicate  unloading. 
Pmax  is  the  maximum  pressure  that  has  been  sustained  by  the  material.  This  parameter 
follows  the  fluid  and  is  used  to  select  the  appropriate  branch  of  the  elastic  curve. 

For  values  of  a  >  1,  the  porous  material  sound  speed  is  defined  as  a  linear  function 
of  porosity: 

,  ,  (a  —  1) 

C=Cs  +  {Ce-Cs)^ - (31) 

(Oo  -  1) 

When  a  =  1,  c  is  computed  from  the  Mie-Grueneisen  relations.  Here  Cs  and  Ce  are  the 
reference  speeds  of  sound  in  the  fully  dense  solid  and  virgin  materials. 


ELASTIC  AND  PLASTIC  a  RELATIONS 

As  shown  in  Figure  1,  the  function  A  has  multiple  branches  and  features  an  unloaded 
(virgin)  value  of  a=ao.  The  appropriate  branch  of  A  is  selected  as  follows: 

1.  When  loaded,  the  virgin  material  compresses  elastically  along  curve  A®,  until  a  pressure 
of  pe  is  reached  (i.e.  point  A). 

2.  In  response  to  further  loading,  the  plastic  curve  Apj  is  followed  until  the  pressure  reaches 
Ps.  Here  a=l,  a  value  which  is  retained  forever. 

3.  The  unloading  path  followed  from  a  point  on  the  plastic  curve,  is  along  the  intersecting 
elastic  curve.  For  example,  point  B  unloads  along  A^j. 

4.  From  a  position  on  an  elastic  curve,  the  loading  path  follows  the  elastic  curve  until  the 
plastic  curve  is  reached.  Further  loading  is  along  the  plastic  curve. 
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The  above  description  of  p-a  material  behavior  indicates  that  a  state  {p,  pmax ,  Pt)  is  on  the 
plastic  curve  if  pt  >  0  and  a  =  Api(p);  otherwise  it  is  on  the  elastic  curve  which  intersects  the 
plastic  curve  at  pmax-  Noting  that  A  is  a  monotone  function  of  p,  it  is  possible  to  simplify 
this  rule:  a  state  is  on  the  plastic  curve  if  p  >pmax  and  on  the  elastic  curve  otherwise. 

Alternatively,  it  is  possible  to  cast  the  test  for  the  plastic  curve  location  in  terms  of 
a:  a  state  is  on  the  plastic  curve  if  a^Omin  and  otherwise  on  the  elastic  curve  intersecting 
o^min-  Here  amin  is  the  minimum  value  of  a  attained  by  the  material  in  the  cell.  Since  it 
is  a  historical  marker,  it  must  be  convected  with  the  fluid  and  it  is  initialized  by  setting 
=  Oe-  This  rule  can  be  expressed: 


^min 

a-  [p^oirmn)  -  |  Aei{p,amin)  otherwise 
Inverting  A,  as  is  shown  in  Figure  2,  and  solving  for  p  yields: 


P  —  P{,^t  ^min) 


Ppl{oi)  if  «  <  OCmin 

Pei{a,amin)  otherwise 


(32) 


(33) 


The  equation  for  the  plastic  curve  is: 

Ap/(P)  =  1  +  (Oe  -  1)  > 

\PS-PeJ 

which  can  also  be  inverted  to  obtain: 

(g-  1)  1^ 

(«e  -  1). 


Pplio^)  =  PS-  {PS  -  Pe) 


(34) 


(35) 


Traditionally,  the  elastic  curve,  Agi,  is  given  in  the  differential  form 

dl^lia)  _o?  (.  1  Y  -  1  I  -c«)(a-  1) 

dp  ~  Kc\  h?)'  Cs{ao-l) 


(36) 


where  Kc  is  the  bulk  modulus  of  the  solid  material.  In  the  present  work,  this 

differential  relation  has  been  replaced  by  a  linear  one: 


Ae/(P)  OLmin) 


(37) 


The  significance  of  this  change  can  be  appreciated  by  constructing  the  P-g  elastic-plastic 
curves  from  Equations  (34)  and  (36)  for  parameter  values  typical  of  sand.  Results  are  shown 
in  Figure  3  and  indicate  that  the  elastic  curves  are  nearly  linear  except  as  g  ^  go  where  the 
elastic  curves  displays  a  nonlinear  behavior,  rising  above  the  plastic  curve,  A^i  . 

.  Inverting  Equation  (37)  and  solving  for  p  yields: 


Pel{oc,amin)  =  Ppl{oimin)  +  ( 

dp 


(38) 
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COMPUTING  p  AND  p 

The  strategy  for  implementing  the  P-a  model  in  DYSMAS  is  to  make  this  it  appear 
similar  to  a  conventional  equation  of  state.  All  the  work  needed  to  support  this  equation  of 
state  is  encapsulated  within  a  single  routine  that  is  called  as  a  function  of  /o,e  and  amin-  This 
is  accomplished  by  equating  Equations  (28)  and  (33)  and  iteratively  solving  the  resulting 
nonlinear  relation 


p{a,p,e)  =  P{o(,amin) 


(39) 


for  a. 

The  solution  of  Equation  (39)  requires  the  selection  of  the  appropriate  branch  of  the 
P{a,  ocjnin)  curve  of  which  there  are  three  possibilities:  elastic,  plastic,  and  the  a=l  extension 
for  compressed  material.  These  three  branches  are  illustrated  in  Figure  2  for  a  =  Also 
shown,  are  three  graphs  for  Equation  (28),  each  associated  with  a  different  solution  branch. 
The  graphs  for  the  p{a^  /^i,  ei),  p{a,  p2, 63),  and  p{a,  pz,  63)  curves  intersects  the  q;=1  line, 
the  plastic  curve  and  the  elastic  curve,  respectively.  The  iterative  procedure  first  determines 
the  appropriate  branch  of  the  P  curve  and  then  calculates  the  point  of  intersection.  The  a=l 
extension  is  selected  if  p(l,/9,  e)  >  ps  or  amin  =  1;  otherwise  the  plastic  branch  applies 
if  p{amin-,  p-,^)  >  Ppi{otmin',oimin)‘  If  both  tests  fail,  the  intersection  with  the  elastic  curve 
Pei{ci,amin)  is  Computed. 

The  computation  of  p  given  amin*  P.  e,  is  a  two  step  process;  compute  a  from  Equation 
(33)  and  then  use  the  Mie-Grueneisen  function  to  determine  the  value  p,  which  satisfies: 

ap  =  ps{ap,e).  (40) 


CHANGES  OF  STATE 


The  isentropic  state  change  is  computed  by  integrating 

P2 


Pi 


For  a  P-q;  material,  a  closed-form  solution  is  not  available  and  this  equation  must  be  solved 
numerically.  To  improve  efficiency,  the  numerical  integration  need  only  be  carried  out  when 
a  >  1.  The  remaining  part  of  the  integration  can  be  computed  from  the  Mie-Grueneisen 
equation  of  state  relation.  Equation  (17). 

Hugoniot  changes  of  state  are  computed  by  solving  the  Hugoniot  condition: 


Pi +P2  = 


2pip2{e2  —  ei) 
{p2  -  pi) 


(42) 
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If  a  I  and  012  are  known.  Equation  (42)  can  be  written  in  terms  of  the  solid  Mie-Gnieneisen 
properties  {ps,Ps)  and  e  can  be  eliminated.  Reintroducing  the  porous  material  properties 
yields: 

,  C\p\p2—p\C2 

P2 —  Pi -r 


(p2-pl)  _  PiPjOl 
2  Topo 


where: 


Cl  =  er{a2P2)  -  er(ai/>i) 


\Pr{ot2p2)  -  p{aipl)\ 


C2  =  {P2  —  pi) 

With  p2  known,  62  follows  from: 


p2pi{<^2  —  «i) 


FoPo 


62  =  6^(025^2)  + 


FoPo 

[q;2P2  -pT{ot2,p2)] 


FoPo 


(43) 


(44) 


(45) 


Since  02  is  often  not  known,  it  may  be  necessary  to  imbed  these  relations  in  an  iterative 
procedure. 


DYSMAS  IMPLEMENTATION 

A  necessary  consideration  in  coupling  the  P-or  equation  of  state  to  DYSMAS  is  com¬ 
patibility  with  the  delta  formulation.  The  solid  state  regime  of  the  P-a  equation  of  state  is 
much  stiffer  than  the  porous  one,  and  consequently,  there  is  a  greater  need  for  accuracy  here. 
Hence  the  solid  density  and  energy  reference  states  are  used  to  normalize  the  density  and 
energy  and  the  following  definitions  apply: 

Ap  —  P  Poo )  ^Ps  —  Ps  Poo 

Ap  =  p  -  poo]  ^Ps  =  Ps  -  Poo  (46) 

^^6  6  ^OOj  ^^^3  ^3  Poo 

Substituting  the  above  into  Equations  (28)  and  (29)  yields  the  following  relations  between 
solid  and  porous  material  delta  properties: 

Aps  -  Ap  +  poo{oi  -  1) 

Aes  =  Ae  (47) 

Aps  =  Ap  +  poo(«-  1)  • 

The  required  DYSMAS  P-a  equation  of  state  functions  are  as  follows: 

•  Function  Fi:  Compute  Api,cf,ai  given  Api,Aei,amm- 

•  Function  F2:  Compute  Ap2,Ae2, 02,02  given  Api,  Aei,  Ap2,OTOm  and  assuming  an 
isentropic  change  of  state  between  Api  and  Ap2- 

•  Function  F3:  Compute  Ap2,Ae2,C2,a  given  Api,  Aei,  Ap2,Omm  and  assuming  a 
change  of  state  along  the  Hugoniot  between  Api  and  Ap2- 

•  Function  F4:  ComputeApi,Cj,Q’i  given  Api,  AeijO^m 

•  Function  F5:  Compute  Ap2,Ae2, 02,02  given  Api,  Aei,  Ap2,ctmm  and  assuming  an 
isentropic  change  of  state  between  Api  and  Ap2. 

•  Function  Fe:  Compute  Ap2,  Ae2,c2,02  given  Api,  Aei,  Ap2, Omm  and  assuming  a 
change  of  state  along  the  Hugoniot  between  Api  and  Ap2. 
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A  detailed  description  of  these  routines  is  provided  in  Appendix  B 

The  interaction  between  the  DYSMAS  scheme  and  the  P-a  equation  of  state  during  a 

computational  step  can  be  summarized  as  follows: 

•  In  phase  1,  calculated  c  from  Equation  (33)  with  a=a/.  Here  a/  is  the  value  of  a 
computed  in  the  previous  cycle. 

•  Execute  the  standard  Euler  algorithm,  phases  1-5.  However,  add  the  convection  of  amin 
using: 

+  Y-V{pamin)  =  0  (48) 

•  Use  the  convected  value  of  Omin  in  phase  6  to  service  equation  of  state  calls  which  then 
return  the  new  value  of  a=ai 

•  At  the  end  of  Phase  6,  for  each  cell  define: 

OCmin  =  min{ai,arnin)  (49) 
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CHAPTER  4 
NUMERICAL  RESULTS 

The  numerical  scheme  outlined  in  the  preceding  section  has  been  applied  to  several  types 
of  cases.  The  first  type  is  designed  to  validate  implementation  of  the  P-a  and  compares 
DYSMAS  to  exact  solutions  and  other  codes  containing  the  P-a  model.  The  second  type 
tests  the  ability  of  the  P-a  equation  of  state  to  simulate  sand.  This  is  accomplished  by 
comparing  calculation  and  experiment. 

VALIDATION  TESTS 

The  implementation  of  the  P-a  model  in  DYSMAS  is  tested  by  comparing  with  an  exact 
Riemann  solution  as  well  as  with  the  CTH^  and  MMG^  codes.  The  Mie-Grueneisen  and  P-a 
parameters  used  in  the  verification  tests  are  listed  in  Table  1. 


QUANTITY 

Verification 

Tests 

SRI  100% 

SRI  95% 

SRI  78%. 

SAMSI 

S 

1.93 

2.17 

2.18 

2.13 

2.2 

Cs 

.148(10^) 

.2(10^) 

.2(10^) 

.2(10^) 

.14(10^) 

Po 

2.070 

2.052 

2.052 

2.012 

1.444 

To 

.880 

.865 

.870 

.846 

.610 

(Xq 

1.052364 

N/A 

1.01988 

1.10854 

1.00958 

Ps 

6.5(10'^) 

N/A 

5.(10^) 

5.(10^) 

.138(10*) 

Ce 

6.0(10^) 

N/A 

6.0(10^) 

6.0(10^) 

1.13(10^) 

Pe 

0. 

N/A 

0. 

0. 

0. 

0:e 

1.052364 

N/A 

1.01988 

1.10854 

1.00958 

TABLE  1.  MIE-GRUENEISEN  AND  P-a  PARAMETERS  FOR  TEST  CASES. 


Riemann  Problem  Solutions 

Solutions  for  fully  saturated  sand  and  partially  saturated  sand  are  considered.  Fully 
saturated  sand  is  represented  by  the  Mie-Grueneisen  Equation  of  state  while  the  P-a  Equation 
of  state  is  used  to  model  the  partially  saturated  sand.  Computed  densities  for  these  two  cases 
are  shown  in  Figures  4  and  5  and  close  agreement  between  computation  and  experiment  is 
achieved.  These  calculations  were  run  using  500  mesh  points.  The  left  and  right  states  for 
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the  fully  saturated  case  were:  pi  =  2.20, e;  =  5(lO^),Pi  =  4.69(10®);  pr  =  2.054, Cr  = 
5(l0®),pr  =  4.42(10^)  Ur  =  0.  The  conditions  for  the  partially  saturated  case  were: 
PI  =  4.23(10®), =  2.20, e/  =  5(l0®);  pr  =  4.42(10^), Pr  =  1.97,6^  =  5(10®).  In  both 
cases,  left  and  right  state  velocities  were  initially  zero. 


1-D  Spherical  bubble  for  Fully  Saturated  Sand 


This  problem  consists  of  a  spherical,  high  pressure  gas  bubble  containing  a  gamma  law 
gas  (7  =  1.3)  imbedded  in  fully  saturated  sand.  The  sand  is  described  with  a  Mie-Grueneisen 
equation  of  state  and  the  computation  was  completed  using  a  CFL  number  of  .8.  The 
computed  spatial  pressure  and  velocity  profiles  are  exhibited  for  the  DYSMAS,  CTH,  and 
MMG  codes  in  Figures  6  and  7  at  five  different  times  during  the  calculation.  The  results 
from  both  codes  are  in  good  agreement  except  inside  the  bubble,  where  the  velocities  differ. 
Temporal  plots  of  the  velocity  at  points  within  the  bubble  exhibit  rapid  oscillation  in  all 
three  of  the  codes,  reflecting  the  repeated  reflection  of  shock  and  expansion  waves  across  the 
bubble,  in  response  to  the  bubble  surface  motion.  To  the  extent  that  the  bubble  surface  does 
not  move  in  precisely  the  same  way  in  each  code,  the  velocities  within  the  bubble  differ. 

1-D  Spherical  bubble  for  Partially  Saturated  Sand 


This  cases  is  similar  to  the  preceding  one  with  the  exception  that  the  sand  is  not  fully 
saturated.  Accordingly,  the  sand  was  modeled  using  a  P-q  equation  of  state  with  ao  taken 
to  be  1.052808.  The  computation  was  completed  using  a  CFL  number  of  .45  and  FCT 
diffusion  and  antidiffusion  values  of  .125. 

The  computed  spatial  pressure,  velocity  and  a  profiles  are  shown  at  four  different  times 
using  the  DYSMAS,  CTH  and  MMG  codes  in  Figures  8,  9  and  10.  As  in  the  fully  saturated 
sand  case,  pressure  agrees  well  everywhere,  while  the  velocities  are  in  good  agreement 
within  the  sand.  The  a  values  are  also  in  close  agreement  everywhere.  Figures  11  and  12 
provide  the  pressure  and  velocity  time  histories  at  three  different  locations  near  the  explosion 
bubbles.  These  figures  compare  the  DYSMAS,  MMG,  and  CTH  codes  and  close  agreement 
is  obtained. 

Figures  13  and  14  demonstrate  the  influence  of  the  FCT  diffusion  and  antidiffusion 
parameters  (FCTDIF  and  FCTADF)  on  the  spatial  pressure  and  velocity  distribution.  Each  of 
these  figures  exhibits  calculations  with  FCTDIF=FCTADF=.125  and  FCTDIF=FCTADF=.5. 
Changes  in  the  FCT  parameters  influence  the  pressure  near  the  shock  and  near  the  bubble. 
Increasing  FCT  values  widened  the  compute  front  and  decreased  the  pressure  near  the  bubble 
surface. 
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2-D  TNT  Explosion 

The  computational  domain  is  shown  in  Figures  15  and  16  and  features  partially  saturated 
sand  in  the  lower  and  water  in  the  upper  half  of  the  domain.  The  initial  TNT  charge  is  located 
on  the  sand  siuface.  The  computation  was  completed  using  a  CFL  number  of  .45  and  FCT 
parameter  settings  of  .5.  The  calculation  was  viable  for  lower  values  of  the  FCT  parameters; 
however,  the  solutions  exhibited  excessive  oscillations.  The  pressure  distribution  is  shown 
in  Figure  15  and  clearly  illustrates  the  propagation  of  the  explosion  shock  through  the  water 
and  sand.  The  diminished  shock  radius  in  the  sand  along  the  vertical  is  a  consequence  of 
the  lower  sound  speed  in  this  material.  The  sand  shock  shape  curves  outwards  near  the 
water/sand  interface.  This  is  a  consequence  of  the  pressure  imposed  on  the  sand  by  the 
shocked  water  near  this  interface. 

Figure  16  illustrates  the  behavior  of  a  in  the  computational  domain.  The  black  level 
is  that  of  the  virgin  sand  (a  =  1.052808)  while  the  light  level  in  the  water  and  behind  the 
sand  shock  represents  a  value  of  1.  It  is  evident  that  the  passage  of  the  shock  compresses 
the  sand  material  to  a  solid  state. 

Computed  pressure  profiles  are  compared  for  the  DYSMAS  and  MMG  codes  in  Figures 
17,  18  and  19.  The  first  of  these  figures  illustrates  the  pressure  profile  along  the  centerline 
of  the  computational  domain.  The  remaining  two  figures  demonstrate  horizontal  traverses 
through  the  flow  field,  slightly  above  and  slightly  below  the  sand/water  interfaces.  This  three 
paths  are  indicated  by  the  dotted  lines  in  Figure  15. 

APPLICATION  OF  P-a  TO  SAND  TESTS 

The  capability  of  the  P-a  equation  of  state  to  model  sand  is  tested  by  comparing  with 
experiment.  Two  different  sets  of  data  are  used:  SRI  1— D  sand  explosions  and  SAMSI  2— D 
explosive  induced  deformation  of  the  buried  target.  In  each  case,  initial  estimates  of  the 
Mie-Grueneisen  equation  of  state  parameters  were  constructed  from  the  water  and  sand  grain 
value  as  described  in  Appendix  E.  Additional  p-a  parameters  were  based  on  test  results  (if 
available)  and  experimentation. 

SRI  Particle  Velocity  Tests 

The  SRI  tests  consist  of  an  explosion  within  a  sand  filled,  cylindrical  container^.  Particle 
velocities  at  different  distance  from  the  explosion  were  measured  as  a  function  of  time.  The 
resulting  experiment  is  1-D  until  the  shock  arrives  at  the  container  walls.  Measurements 
were  made  using  fully  saturated,  95%  saturated  and  78%  saturated  sand. 

The  Mie-Grueneisen  parameters  (Cj,  S,  and  To)  for  the  fully  saturated  sand-water 
mixture  were  adjusted  to  best  match  experiment.  The  initial  density  was  available  from 
the  experiment.  For  partially  saturated  sand,  additional  parameters  need  to  be  specified:  ao, 
ps,  and  pe-  OLo  was  computed  from  the  measured  dry,  fully  saturated  and  partially  saturated 
sand  density,  pe  was  taken  as  zero  and  P*  was  adjusted  to  best  match  experiment.  The 
selected  set  of  parameters  is  given  in  Table  1. 
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The  measured  results  for  the  SRI  test  are  shown  in  Figures  20,  21  and  22  along  with 
the  computational  results  for  each  test.  An  examination  of  these  figures  indicates  that  the 
computed  velocity  immediately  following  the  arrival  of  the  explosion  shock  are  in  reasonable 
agreement  with  experiment.  However,  experimental  results  diminish  faster  than  the  calculated 
results. 

SAMSI  Target  Deformation  Tests  Conducted  At  WES 

In  each  of  these  three  tests^,  a  mine-like  target  was  buried  in  a  sand-like  soil  and  subjected 
to  an  explosive  shock.  The  experimental  arrangement  for  this  test  is  depicted  in  Figure  23a 
and  features  8  lbs  of  pentolite  and  an  Sin  in  diameter  aluminum  target.  Pressures  were 
monitored  at  several  locations  in  the  sand  and  the  final  deformation  of  the  mine  face  was 
recorded. 

The  first  SAMSI  test  has  been  simulated  with  a  coupled  Eulerian/Lagrangian  DYSMAS 
calculation.  The  setup  for  this  calculation,  including  the  grids  for  the  Euler  and  Lagrangian 
modules,  are  depicted  in  Figure  23b.  The  calculation  was  performed  using  2-D,  cylindrical 
coordinates,  with  the  coordinate  axis  orientated  along  the  center  of  the  explosive  charge  and 
the  mine.  The  water  above  the  sand  and  the  in-situ  soil,  located  3  ft  below  the  mine  were 
ignored. 

The  “sand”  in  this  test  consisted  of  a  grout  with  a  nominally  .5%  to  1%  air  content. 
Density  was  the  only  piece  of  information  available  to  characterize  the  sand.  The  remainder 
of  the  P-a  equation  of  state  parameters  were  determined  empirically  by  comparing  calculation 
and  sand  pressure  measurements.  Table  1  provides  the  selected  set  of  parameter  values. 

The  evolution  of  the  flow  field  is  shown  in  Figures  24  to  27.  Here  the  passage  of  the 
shock  through  the  sand  is  clearly  visible  as  is  the  interaction  between  the  shock  and  the 
deforming  target  surface.  The  weak  shock  visible  in  Figure  26  above  the  cylinder  arises 
from  the  initial  impact  shock  transmitted  through  the  target.  A  comparison  of  the  measured 
and  the  calculated  pressures  in  the  sand  is  provided  in  Figure  28.  The  calculations  are  in 
reasonable  agreement  with  experiment,  although  the  computed  pressure  level  is  marginally 
higher  than  experiment  for  the  gauge  near  to  the  explosion. 

Figure  29  presents  the  displacement  history  at  the  target  surface  centerline.  Also  included 
in  this  figure  is  the  target  translation  and  a  second  curve  which  is  the  difference  between 
the  displacement  and  the  translation.  This  final  curve  represents  the  centerline  deformation 
The  final  computed  target-surface  deformation  profile  is  shown  in  Figure  30  and  is  in  good 
agreement  with  measured  values,  which  are  also  given. 
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CHAPTER  5 

SUMMARY  AND  CONCLUSIONS 

The  Mie-Grueneisen  and  P-a  equations  of  state  have  been  incorporated  into  the  DYSMAS 
code.  The  form  of  these  equations  of  state  follows  from  Reference  3;  however,  the  manner 
of  implementation  has  been  altered  to  provided  compatibility  with  the  DYSMAS  code.  In 
particular,  these  equations  were  cast  in  the  DYSMAS  delta  form  and  two  auxiliary  functions 
were  constructed  to  compute  isentropic  and  Hugoniot  jumps.  In  addition,  a  second  form  for 
these  equations  of  state  and  auxiliary  functions  were  developed  which  were  a  function  of 
(p,e)  rather  than  the  usual  (p,e).  The  inverted  equation  of  state  form  and  auxiliary  function 
were  implemented  numerically. 

The  P  -a  equation  of  state  routine  was  constructed  to  be  similar  to  other  equations  of 
state  and  differs  only  in  the  need  to  specify  the  minimum  a  .  This  parameter  represents 
the  minimum  value  of  a  attained  by  the  material  at  the  point  of  interest.  Thus,  the  form 
for  the  equation  of  state  is  p(a:min./o,e)  or  /jComin.P.e).  A  Newton  iteration  was  performed  to 
determine  the  pressure  or  density  associated  with  the  calling  arguments. 

The  implementations  of  the  Mie-Grueneisen  and  P-a  equations  of  states  have  been 
numerically  verified  by  comparing  with  exact  solutions  and  the  results  from  other  codes. 
The  exact  solutions  cases  were  1-D  planar  Riemann  problems,  while  the  code  comparison 
cases  consisted  of  a  spherical  explosion  (explosive  surround  by  sand)  and  a  2-D  simulation 
of  an  imderwater  explosive  sitting  on  a  sand  bottom.  In  aU  cases,  the  DYSMAS  results 
compared  closely  to  the  appropriate  benchmark. 

The  veracity  of  P-a  as  a  sand  model  was  tested  by  comparing  with  two  experiments. 
The  first  was  the  SRI  explosion  in  a  canister  test  which  was  simulated  using  a  1-D  spherical 
calculation.  The  second  test  featured  a  target  buried  in  sand  that  was  subjected  to  an 
explosion.  In  the  latter  case,  the  Euler  code  was  coupled  to  the  Lagrangian  code,  and 
target  deformation  as  well  as  sand  pressure  were  compared  to  experiment.  In  both  cases  the 
predictions  and  experiment  were  in  reasonable  agreement  with  experiment. 

»  The  numerical  results  presented  in  this  paper  establish  the  validity  of  the  P-a  implemen¬ 
tation  in  the  DYSMAS  code.  In  addition,  agreement  between  calculation  and  experiment 
demonstrated  on  the  two  test  cases  is  encouraging.  However,  additional  comparisons  with 
experiment  are  necessary  to  more  fully  ascertain  the  merits  of  the  P-a  equation  of  state  as 
a  sand  model. 
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GLOSSARY 


symbol 

ao...a5. 

A 

c 

Ce 

Cs 

e 

Co 

er 

p 

Pe 

Pmax 

Po 

Pr 

Ps 

PS 

Pt 

P 

S 

S 

Y 

a 

Oe 

ai 

O-o 

^min 

Ap,Ap,Ae 


Definition 

see  Equation  (6) 

a  function,  Equation  (32) 

Sound  speed 

Reference  sound  speed  in  porous  material 
Reference  sound  speed  in  solid  (non-porous)  material 
Energy 

Mie-Gruneisen  equation  of  state  reference  state  energy 
Mie-Gruneisen  equation  of  state  reference  energy 
Pressure 

Pressure  plastic  behavior  of  P-a  material  is  initiated 

Maximum  pressme  experienced  experienced  by  material 

Mie-Gruneisen  equation  of  state  reference  state  pressure 

Mie-Gruneisen  equation  of  state  reference  pressure 

Solid  material  pressure 

Pressure  at  which  all  pores  are  crushed 

Material  derivative  of  pressme 

Inverse  of  A  function.  Equation  (33) 

Entropy 

Coefficient  relating  shock  (uo)  and  particle  velocity  (Up):  Ug  =  Cq-V  Sup 
See  Equation  (5) 

Degree  of  porosity,  see  Equation  (29) 

Maximum  a  where  Api  is  valid 
Local  a  at  a  particular  cell 
Initial  a  at  p=0 

Minimum  value  of  a  achieved  by  material 
See  Equation  (21) 
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Apr,Aer 

To 

P 

po 

Ps 


See  Equation  (21) 

See  Equation  (46) 

Gruneisen  parameter  at  po 
Equation  (4) 

Density 

Mie-Gruneisen  equation  of  state  reference  state  density 
Solid  material  density 


Subscripts 

1,  2  Initial  and  final  equation  of  state  conditions 

el  Elastic  curve 

pi  Plastic  curve 

s  Isentropic  or  solid 
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FIGURE  3.  NONLINEAR  P-a  ELASTIC-PLASTIC  CURVES  FOR  A  SAND  LIKE  MATERIAL. 
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nCURE  4.  COMPUTED  AND  EXACT  DENSITY  FOR  THE  MIE-GRUENEISEN  RIEMANN  PROBLEM 
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FIGURE  5.  COMPUTED  AND  EXACT  DENSITY  FOR  THE  P-a  RIEMANN  PROBLEM 
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FIGURE  6.  PRESSURE  RESPONSE  OF  FULLY  SATURATED  SAND  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE 
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HGURE  7.  VELOCITY  RESPONSE  OF  FULLY  SATURATED  SAND  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE 
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HGURE  8.  PRESSURE  RESPONSE  OF  PARTIALLY  SATURATED  SAND  (ao=L052808)  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE 
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FIGURE  9.  VELOCITY  RESPONSE  OF  PARTIALLY  SATURATED  SAND  (q„=L052808)  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE 
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HGURE  10.  COMPUTED  a  FOR  PARTIALLY  SATURATED  SAND  (ao=l. 052808)  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE 
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FIGURE  11.  PRESSURE  HISTORIES  AT  THREE  POINTS  FOR  PARTIALLY  SATURATED 
SAND  (ao=  1.052808)  IN  RESPONSE  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE 


FIGURE  12.  VELOCITY  HISTORIES  AT  THREE  POINTS  FOR  PARTIALLY  SATURATED 
SAND  (a„=L052808)  IN  RESPONSE  TO  A  GAMMA  LAW  EXPLOSION  BUBBLE 
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FIGURE  13.  INFLUENCE  OF  FCT  PARAMETERS  ON  THE  COMPUTED  PRESSURE  FOR  A  GAMMA  LAW  EXPLOSION  BUBBLE. 
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FIGURE  14.  INFLUENCE  OF  FCT  PARAMETERS  ON  THE  COMPUTED  VELOCITY  FOR  A  GAMMA  LAW  EXPLOSION  BUBBLE. 
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FIGURE  17.  PRESSURES  ALONG  THE  CENTERLINE  OF  THE  EXPLOSION  OF  FIGURES  15  AND  16. 


le+10 


NSWCDDyTR-95/107 


(3¥¥UI0/P) 


d 


FIGURE  18.  PRESSURES  ALONG  THE  z=-550  LINE  OF  THE  EXPLOSION  OF  RGURES  15  AND  16. 
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FIGURE  19.  PRESSURES  ALONG  THE  z=-450  LINE  OF  THE  EXPLOSION  OF  FIGURES  15  AND  16. 
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nOURE  20.  COMPUTED  AND  MEASURED  VELOCITffiS  FOR  THE  100%  SATURATED  SRI  TEST. 
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FIGURE  21. .  COMPUTED  AND  MEASURED  VELOCITIES  FOR  THE  95%  SATURATED  SRI  TEST. 
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HGURE  24.  COMPUTED  PRESSURES  FOR  THE  SAMSI  TEST  AT  734  //SECS. 


nOURE  25.  COMPUTED  PRESSURES  FOR  THE  SAMSI  TEST  AT  789  //SECS. 


41 


NSWCDD/TR-95/107 


FIGURE  26.  COMPUTED  PRESSURES  FOR  THE  SAMSI  TEST  AT  844  fxSECS. 


FIGURE  27.  COMPUTED  PRESSURES  FOR  THE  SAMSI  TEST  AT  953  /zSECS. 
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Fig.  28  SAMSI  Calibration  Shot 


FIGURE  28.  COMPUTED  AND  MEASURED  SAND  PRESSURES  FOR  THE  SAMSI  TEST. 
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FIGURE  29.  COMPUTED  CENTERLINE  DISPLACEMENT  HISTORY  FOR  THE  SAMSI  TEST. 
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HGURE  30.  COMPUTED  AND  MEASURED  TARGET  DEFORMATION  t^OR  THE  SAMSI  TEST. 
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APPENDIX  A 

MIE-GRUENEISEN  EQUATION  OF  STATE  ROUTINES 
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This  appendix  provides  a  description  of  the  MieGrueneisen  equation  of  state  subroutines. 

1.  Subroutine  STDMIE:  /i(A/)i,Aei)  (Api,Ci,  Ap^,  AcrJ 

a.  Function:  Compute  Api,Cj  given  Api,Aei. 

b.  Description:  Evaluate  appropriate  equations. 

c.  Algorithm:  Solve  Eq.  (22)  for  Ap  and  Eq.  (23)  for  c^.  Evaluating  these  equations 
requires  use  of  Eqs.  (2)-(6),  (11)  and  also  produces  Apn,  Ae^. 

2.  Subroutine  ISEMIE:  /2(Api,  Aei,  Ap2)  (Ap2,  Ae2,C2) 

a.  Function:  Compute  Ap2,Ae2,C2  given  Api,Aei,p2  and  assuming  an  isentropic 
change  of  state  between  Api  and  Ap2- 

b.  Description:  Evaluate  Equation  (25  ) 

c.  Algorithm: 

/i(Api,Aei)  (Api,cf, Apn,  Ae^i) 

/i(Ap2,Aei)  -+  (Ap«,c<,Apr2,  Acrj) 

compute  Acz  from  Eq.  (16)  and  solve  Equation  (22)  for  Ap2. 

Re-evaluate  c:  =  cf  -f  roPo(Ap2  -  Apt)/p| 

3.  Subroutine  HUGMIE:  /3(Api,  Aei,  Ap2)  ->  (Ap2,Ae2,c|) 

a.  Function:  Compute  Ap2,Ae2,C2  given  Api,Aei,Ap2  and  assuming  a  change  of 
state  along  the  Hugoniot  Api  and  Ap2. 

b.  Description:  Evaluate  Equation  (26). 

c.  Algorithm: 

/i(Api,Aei)  *  (Api , Cj,  Aprj ,  Acfj) 

/i(Ap2,Aei)  (Api,Q,Apr2,  Ae^J 
Compute  5p  from  Eq.  (26). 

Ap2  =  Api  -f-  6p 

Solve  Eq.  (22)  for  Ae2. 

Re-evaluate  c:  +  Topo{Ap2  -  Apt)/p| 

4.  Subroutine  ISTMIE:  /4(Api,Aei)  (Api,c^) 

a.  Function:  Compute  Api,Cj  given  Api,  Aei 

b.  Description:  Iterate  for  the  value  of  p  which  statistics  Equation  (22). 

c.  Algorithm: 

ApO  =  0 

Iterate  from  n=l  until  convergence  or  n  >  maximum  number  of  iteration. 
/i(Ap”-\Ae)  (Ap”,c2") 

compute  from  Eq.  (12). 

Compute  8p^  =  (Ap2  -  Ap”)/(|^)^ 
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6 <  10“^  iteration  is  complete.  Otherwise: 

Ap”  =  +  min\max{8 ,  “-2),  .2]  (1) 


Set:  Api  =  A/)";  c  = 


5.  Subroutine  IISMIE:  /5(Api,  Aei,  Ap2)  (Ap2,Ae2,C2) 

a.  Function:  Compute  p2,e2,cl  given  pi,ei,p2  and  assuming  an  isentropic  change  of 
state  from  Api  and  Ap2- 

b.  Description:  Apply  ISTMIE  to  determine  initial  density.  Iterate  for  a  final  value  of 
rho  which  is  on  the  same  isentrop  as  the  initial  conditions  and  has  a  pressure  value 
equal  to  p2 

c.  Algorithm: 

fi{Api,Aei)  (Api,cf) 

=  (Ap2  -  Api)/c^;  Ae2  =  Aei;Ap^  =  Api 
Iterate  from  n=l  until  convergence  or  n  >  maximum  number  of  iteration. 

A/o”  =  A/o”"^  +  min[max{6p''’’,—.2),  .2] 

/i(Ap",Ae”-i)  (ApuAer,,cj) 

Solve  Eq.  (25)  for  Ae" 

Solve  Eq.  (22)  for  Ap" 

Re-evaluate  Ct:  +  ToPoiAp^  —  Apt)/p2 

8p^  =  {Ap^-Ap^)lc\ 

If  8p^  <  10“"^  iteration  is  complete. 

Ap2  =  Ap^-,  Ae2  =  Ae” 

6.  Subroutine  IHUMIE:  /6(Api,  Aei,  Ap2)  — »•  {Ap2,Ae2,C2) 

a.  Function:  Compute  p2,e2,cl  given  pi,ei,p2  and  assuming  a  change  of  state  along 
the  Hugoniot  from  Api  and  Ap2. 

b.  Description:  Apple  ISTMIE  to  determine  the  initial  density.  Iteratre  for  a  final  value 
of  denity  which  as  along  Hugoniot  and  produces  the  pressure  p2 

c.  Algorithm: 

/3(Api,Aei)  (Api,c^) 

Iterate  from  n  =  1  until  convergence  or  n  >  maximum  number  of  iterations.  Start 
iteration  with:  8p^  =  8p/cl]  Ae”  =  Aei;  Ap”  =  Api 

Ap”  =  Ap”“^  -f  mm[maa;(6p”,  — .2),  .2] 

Solve  Eq.  (17)  for  Ae^  and  set  Ae”  =  Ae^ 

/i(Api,Aei)  (Ap”,Aer2cf) 

8p  =  {Ap2  -  Ap”)/cf 
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If  <  10  ®,  iteration  is  complete.  Otherwise  : 


8p-  =  -6pl 
8p^  =  -8p 


(  8p-8p-^  \ 

Y  A/?”  --  j 


if  n  >  I 


n  =  l 


A/?2  =  A/)i  +  ;  Ae2  =  Ae” 


A-5 


NSWCDD/TR-95/107 


APPENDIX  B 

P-a  EQUATION  OF  STATE  ROUTINES 
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This  appendix  provides  a  description  of  the  P-a  equation  of  state  subroutines.  In  many 
cases  the  Mie-Grueneisen  Equation  of  State  are  called 

1.  Subroutine  STDP_A:  Fi{Api,Aei,  amin)  (Api,  cf,  Apr, ,  Ae^,  ai) 

a.  Function:  Compute  Api,cf,a:i  given  A/3i,Aei,aTOm- 

b.  Description:  Solves  Equation  (39)  Iteratively 

c.  Algorithm: 

Check  for  the  solid  case.  If  p(l,/o,  e)  >  ps  or  amin  =  1  then  material  is 
compressed  to  a  solid  and  p  =  p{l,p,e)  and  a  =  1 

Otherwise,  check  plastic  case.  If  p{amin,P,e)  >  Pj,i{ot.min,oimin),  use  Newton 
iteration  to  solvep(a, p,  e)  =  Ppi{p,amin)- 

Else  the  case  is  elastic.  Apply  a  Newton  iteration  to  solve  p(a,p,e)  = 
PelijPt  ^min)’ 


2.  Subroutine  ISEP_A:  F2(Api,  Aei,  Ap2,  amin)  (Ap2,  Ae2  c^,  02) 

a.  Function:  Compute  Ap2,Ae2, 02,0:2  given  Api,Ae\,p2,amin  and  assuming  an 
isentropic  change  of  state. 

b.  Description:  Evaluate  Equation  (42)  numerically  using  a  2  point  Runge  method.  Use 
Equation  (25)  for  parts  of  the  change  of  state  path  where  a=l. 

c.  Algorithm: 

Fl(Api,  Aei,Omin)  — >  (Api,Cj,  Apn,  Acri) 

If  Omin^l  CALL  ISEMIE 

Else  to  numerically  ntegrate  Equation  (42). 

apply  a  2  point  Runge-Kutta  method, 

If  0=1  after  any  step  CALL  ISEMIE  to  complete  integration. 

3.  Subroutine  HUGP_A:F3(Api,  Aej,  Ap2,  amin)  (Ap2,  Ae2  c|,  02) 

a.  Function:  Compute  Ap2,Ae2, 02,02  given  Api,  Aei,  Ap2,a„i,„  and  assuming  a 
change  of  state  along  the  Hugoniot. 

b.  Description:  Use  iteration  to  determine  an  02  which  satisfies  Equation  (44). 

c.  Algorithm: 

If  amin  =  1  Solve  for  final  state  with  Equations  (44)-(46)  using  oli  =  a2  =  1. 
Perform  Newton  iteration  ono2  for  final  state  using  Equations  (44)-(46). 
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4.  Subroutine  ISTP_A:  jF'4(Api,  Aei  ?  ^min  )  {Api,cl,ai) 

a.  Function:  Compute  A/9i,Ci,ai  given  Api,Aei,a^,n 

b.  Description:  Compute  a  using  Equation  (33)  and  the  find  the  density  by  calling 
subroutine  ISTMIE 

c.  Algorithm: 

li  p  >  Ps  or  a^nin  =  1  call  ISTMIE  to  determine  pi  and  c^. 

Otherwise: 

Compute  OL  from  Equation  (33)  and  Aps  =  aAp  +  (a  —  l)poo 

/4(Aps,Aei)  ^  {Apa,c]) 

Compute  Ap  =  ^  +  /)oo(J  -  l) 


5.  Subroutine  nSP_A:  F5(Api ,  Aei,  Ap2,ocmin)  {Ap2,  Ae2,  c^,  02) 

a.  Function:  Compute  p2,e2,C2,a2  given  pi, 61,^2,  Q;mm  and  assuming  an  isentropic 
change  of  state. 

b.  Description:  Numerically  integrate  Equation  (42)  from  the  initial  state  until  target 
pressure  is  reached. 

c.  Algorithm: 

•  F4(Api,Aei  5  ^min  )  (Api,cf) 

•  Integrate  Equation  (42)  until  p^2 

Apply  a  2  step  Runge  Kutta  method. 

Fi(Ap”,  Ae”“^a,„m)  (Ap”,cf) 

•  Interpolate  for: 


A/92  =  Ap”-i)(Ap”-Ap” 
Ae2  =  Ae^-d( 


where  : 


Ap2  —  Ap”  ^  \ 
Ap”  —  Ap”“^ ) 


■) 


(3) 


•  Compute  final  properties:i^i(Ap2,Ae2,  Omin  )  {Ap2,cl) 
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6.  Subroutine  IHUP_A:  Fe{Api,Aei,Ap2,amin)  {Ap2,Ae2,cl,a2) 

a.  Function:  Compute  p2,  €2^02,02  given  pi,ei,p2,amm  and  assuming  a  change  of 
state  along  the  Hugoniot. 

b.  Description:  Compute  ai,  02  using  pi ,  p2  and  Equation  (33).  Use  a  Newton  iteration 
to  solve  Equation  (44). 

c.  Algorithm: 

•  If  amin=l.  material  is  fully  compacted.  Compute  p2, 63,  c|  from  feiApi ,  Aei ,  Ap2) 
(Ap2,  Ae2) 

•  Otherwise 

Compute  ai,  0-2  using  pi,  p2  and  Equation  (33). 

Iterate  for  value  of  p2  which  satisfies  Equation  (33) 
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APPENDIX  C 
USER  MANUAL  PAGES 
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Name:  MP_ALP 

Indicates  the  number  of  P-a  materials  included  in  the  calcualtion. 

Default:  MP_ALP=  0 

Restrictions:  MP_ALP>0 

Recommendation:  none 


Remarks: 


Up  to  9  P-a  materials  can  be  included  in  the  calculation.  The  increased  storage  per  cell  is 
2*MP_ALP.  The  reference  and  initial  property  for  a  P-a  material  is  described  in  the  same 
manner  as  other  materials.  Materials  91-95  are  reserved  for  P-a  equations  of  state  while 
96-99  are  for  Mie-Grueneisen  Equation  of  state  materials 
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Codeword:  EQSTPAR  (7/9) 


Mie-Gmeneiseti  Equation  of  State  (96^M<99) 

P  =  Pr{p)-\-^oPo{^- er{p))  (4) 

See  NSWC/TR-95/107  for  definition  of  constants  in  equation. 

Input  parameters 


pos 

variable 

internal  name 

accepted  if 

3 

S 

ESCAPA 

always 

4 

ESCAPE 

always 

5 

Po 

ESD 

always 

6 

Po 

ESALPH 

always 

7 

^0 

ESBETA 

always 

8 

To 

ESGAMM 

always 

9 

Pcav 

ESPCAV 

!=0 

10 

A/>nriax 

DRHOMX 

!=0 

11 

DRHOMN 

!=0 

Remarks: 


The  EQSTPAR  section  should  preceed  the  MATERIAL  section.  To  avoid  confusion 
with  existing  material  definitions,  use  NM=99  for  Mie-Grueneisen  materials  defined  using 
EQSTPAR  section 
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Codeword:  EQSTPAR  (8/9) 


P-a  Equation  of  State  (90<NM^5) 

See  NSWC/TR-95/107  for  a  definition  of  the  P-a  Equation  of  State. 


Input 

parameters 

pos 

variable 

internal  name 

accepted  if 

3 

S 

ESCAPA 

always 

4 

c2 

ESCAPE 

always 

5 

Po 

ESD 

always 

6 

po 

ESALPH 

always 

7 

ESBETA 

always 

8 

To 

ESGAMM 

always 

9 

Q(o 

ESDELT 

always 

10 

Aps 

ESES 

always 

11 

Ce 

ESESP 

always 

12 

Ape 

ESIOTA 

always 

13 

Pcav 

ESPCAV 

!=0 

14 

^Pmax 

DRHOMX 

!=0 

15 

^Pmm 

DRHOMN 

!=0 

Remarks: 


The  EQSTPAR  section  should  preceed  the  MATERIAL  section.  To  avoid  confusion 
with  existing  material  definitions,  use  NM=95  for  Mie-Grueneisen  materials  defined 
using  EQSTPAR  section 
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APPENDIX  D 

STORAGE  ARRYAS  FOR  EQUATIONS  OF  STATE  PARAMETERS 


D-l/D-2 
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The  Mie-Grueneisen  and  P-a  equations  of  state  make  use  of  a  number  of  amterial 
constants.  The  following  table  provides  a  list  of  these  variables  using  the  symbols  defined 
in  the  Glossary. 


Internal  Name 

P-a  constants 

Mie-Gruneisen  constant 

ESCAPA 

S 

S 

ESCAPE 

c2 

‘'S 

ESCAPC 

ai 

at 

ESCAPD 

a2 

a2 

ESA 

as 

as 

ESB 

a4 

a4 

ESC 

as 

as 

ESD 

Po 

Po 

ESALPH 

po 

ESBETA 

eo 

Co 

ESGAMM 

To 

To 

ESDELT 

Oo 

ESES 

Aps 

ESESP 

Ce 

ESETA 

Ape 

ESIOTA 

0:e 

ESPCAV 

Pcav 

DRHOMX 

^/^max 

DRHOMN 
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APPENDIX  E 

COMPUTING  THE  P-a  PARAMETERS  FOR  SAND 


E-l/E-2 
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Sand  consists  of  sand  grain,  water,  and  air.  The  air  and  water  form  a  mixture  which  fills 
the  pore  spaces  between  the  sand  grains.  The  solid  component  of  the  P-a  model  represents 
the  sand  grain  and  water  mixture,  while  the  void  simulates  the  air.  To  apply  the  P-a  model 
to  sand,  the  Mie-Grueneisen  equation  of  state  parameters  for  this  solid  component  must  be 
constructed  from  the  known  water  and  sand  grain  properties. 

The  information  required  to  construct  the  Mie-Grueneisen  constants  for  the  grain-water 
mixture  are: 

1.  The  dry  sand  density,  pd 

2.  The  saturation,  s. 

3.  The  following  material  properties  for  water  and  sand  grain: 

a.  Density,  p. 

b.  Sound  speed,  c. 

c.  The  MieGrueneisen  Constants,  Fq. 

d.  The  shock  velocity  verses  particle  velocity  slope,  S. 

e.  Specific  heat,  Cy  (optional  for  many  applications). 


The  porosity,  n,  which  is  the  volume  fi"action  of  the  void,  is  given  by: 

(Pg  -  Pd) 


n 


Pg 


The  density  of  the  compacted  sand  with  all  of  the  air  pore  volume  eliminated  is: 

_  (1  -  n)pg  4-  snpy, 

Pw+g  —  j:.  r  , 

1  —  (1  —  s)n 

while  the  mass  fractions  of  the  grain  and  water  components  of  this  mixtures  are 


A, I,  — - 


STipxv 


and 


(1  -  n)pg  -I-  snpv, 
\  =  (1 

^  (1  —  n)pg  -f  snpy,  ' 


E-  1 


E-  2 


E-  3 


E-  4 


Furthermore,  from  Equations  (E-2)  to  (E-4)  it  follows  that 

~  ^g^g  ?  E-  5 

where  v  is  the  specific  volume. 

The  specific  heat  for  the  grain-water  mixture  follows  directly  from  a  mass  average: 

Cvg^.u,  —  ^gCvg  +  ^wCvu,  E-  6 

The  speed  of  sound  for  the  grain- water  mixture  arises  from  the  definition: 


ov 


( 

\^p)s 


E-  7 


E-3 
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Assuming  that  the  water  and  grain  components  are  always  in  pressure  equilibrium,  the 
denominator  of  Equations  (E-5)  can  be  evaluated  by  differentiating  Equation  (E-7),  which 
yields: 

dv\  Vg-^-yj  ,  ^  dVp^  ^  f  dVy] 


dp/ ,  c“ 


--A 

-  dp 


-A. 


dp 


_  ,  ,  A  ^ 

—  ^Cf  •>  "r  0 


'3  2 
‘'<7 


E-  8 


Solving  for  the  grain-water  sound  speed  yields: 


^g+w 


^g+w 


The  grain-water  F  also  follows  directly  from  the  definition: 

-o: 

Noting  that  Cg+w  =  \eg  +  and  differetiating  with  respect  to  p  yields: 

“^g-i-w  _  1 


E-  9 


E-  10 


Tg+to  — 


^  _L  ^  ,  A,,,  A 

”  ^  ■+■  v^pn,) 


E-  11 


The  final  property  of  interest  is  S,  the  slope  of  the  Up  verses  Us  relation.  This  quantity 
is  treated  in  a  manner  analogous  to  sound  speed.  Thus, 


C2  _ 


1 


pl+y’i'^p  +  'ik;) 


E-  12 


E-4 
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